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Abstract. We present a new general relativistic hydrodynamics code specifically 
designed to study magneto-rotational, relativistic, stellar core collapse. The code 
is an extension of an existing (and thoroughly tested) hydrodynamics code, which 
has been applied in the recent past to study relativistic rotational core collapse. 
It is based on the conformally-flat approximation of Einstein's field equations 
and conservative formulations for the magneto-hydrodynamics equations. As a 
first step towards magneto-rotational core collapse simulations the code assumes 
a passive (test) magnetic field. The paper is focused on the description of the 
technical details of the numerical implementation, with emphasis on the magnetic 
field module. A number of code tests are presented and discussed, along with a 
representative core collapse simulation. 
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1. Introduction 

During the last decades many efforts have been taken to numerically simulate the 
collapse of stellar cores (see e.g. [TJ [H |3] and references therein). A myriad of 
difficulties arise for its modelling as many physical effects are involved in a process 
whose dynamics is highly nonlinear. Analytic approaches are of limited use and 
numerical simulations are needed to solve the coupled system of equations describing 
the fluid motion, the spacetime evolution, the dynamics of the magnetic fields, and 
the neutrino transport. Such simulations are driving progress in the field despite the 
limited knowledge on relevant physical issues such as realistic precollapse stellar models 
(including rotation and magnetic field strength and distribution) or microphysical, 
finite-temperature equation of state (EOS), as well as numerical limitations involved in 
the challenging task of accounting for Boltzmann neutrino transport, multidimensional 
hydrodynamics, and relativistic gravity. 

Only very recently the first multidimensional simulations of relativistic rotational 
core collapse have become possible, thanks to the use of conservative formulations 
of the hydrodynamics equations and advanced numerical methodology, as well 
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as long-term stable formulations of Einstein's equations (or accurate enough 
approximations) [H[S] (see also [5]). On the other hand, the incorporation of magnetic 
fields and the MHD equations in numerical codes to further improve the realism of 
such simulations in general relativity, is currently an emerging field where significant 
progress is bound to be soon achieved (see [TJ [8] ) . 

Neutron stars have intense magnetic fields (~ fO 12 — fO 13 G) or even larger 
at birth. The presence of such magnetic fields renders magneto-rotational core 
collapse simulations mandatory. In recent years, an increasing number of authors 
have performed axisymmetric simulations (within the so-called ideal MHD limit) 
employing Newtonian treatments of the magneto-hydrodynamics, the gravity, and 
of the microphysics (see [9] and references therein) . The weakest point of all existing 
simulations to date is the fact that the strength and distribution of the initial magnetic 
field in the core are basically unknown. The available simulations show that for weak 
initial fields (< 10 11 G, which is the most relevant case, astrophysics-wise) there are 
no major differences in the collapse dynamics nor in the resulting gravitational wave 
signal, when comparing with purely hydrodynamical simulations. However, strong 
initial fields (> 10 11 G) manage to slow down the core efficiently (leading even 
to retrograde rotation in the proto-neutron star (PNS)) which causes qualitatively 
different dynamics and gravitational wave signals. 

This paper presents a general relativistic magneto-hydrodynamics (GRMHD) 
code designed to study magneto-rotational, relativistic, stellar core collapse. Our 
GRMHD code is an extension of the hydrodynamics code developed by [TOl Hj, 
where magnetic fields are included following the approach laid out in [ITj . Einstein's 
equations are formulated using the conformally flat condition (CFC), which has proved 
very accurate for studying rotational core collapse [T2]. Leaving aside radiation 
transport, whose effects on core collapse dynamics and gravitational radiation are 
dramatic and its numerical modelling is a challenge in itself (see [3] and references 
therein), we adopt here as a first step towards simulating realistic magneto- rotational 
core collapse the passive magnetic field approximation, which is justifiable for the 
magnetic field values present in this scenario. The paper focuses on describing the 
details of the numerical schemes we use and the tests performed with the code. A 
representative magneto-rotational core collapse simulation is also presented and briefly 
discussed. A parameter-space survey of the magneto-rotational core collapse scenario 
will be presented elsewhere. We note that during the development of this work we 
have been made aware of a similar study reported in |f 3] . 

2. Mathematical and physical framework 

We adopt the 3 + f formalism to foliate the spacetime into spacelike hypersurfaces. In 
this formulation the metric reads ds 2 = g^dx^dx v = —a 2 dt 2 + -fij(dx l + fi l dt){dx : ' + 
ftdt), where a is the lapse function, (3 l the shift vector, and jij the spatial three-metric 
induced in each hypersurface. From the energy-momentum tensor T^ v it is possible 
to build the following quantities: E = n^n^T^ = a 2 T 00 , S l = - _Lf n u T flv = 
— i(Toi — Tijft), and Sij =_Lf_LJ — TJy, using the projection operator _L£J and 
the unit four-vector n M normal to the hypersurface. 
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2.1. The CFC approximation 

To solve Einstein's field equations we use the CFC approximation [14] . In this 
approach, the three-metric in the ADM gauge 7^ = (f^jij + hj^ , is assumed to 
be conformally flat, 7^ = 4> 4 %j. In these expressions <j> is the conformal factor, hj^ 
is transverse and traceless, and 7^ is the flat three-metric. Note that the ADM 
gauge choice implies the maximal slicing condition in which the trace of the extrinsic 
curvature vanishes. Under the CFC approximation, Einstein's field equations can be 
written as a system of five elliptic equations 



S = 7 U Sij, and K l i is the extrinsic curvature. 
2.2. General relativistic magneto-hydrodynamics 

For a perfect fluid endowed with an electromagnetic field the energy-momentum tensor 
is the sum of a fluid part and of a electromagnetic field part 



where p is the rest-mass density, h = 1 + e + P/p is the relativistic enthalpy, e is 
the specific internal energy, P is the pressure, and ii M is the four-velocity of the fluid. 
The electromagnetic tensor F» v = U^E y - U U E^ - e^^UxBs, and its dual *F> iU 
can be expressed in terms of the electric field — F^JJv and the magnetic field 
B^ = *F^ V U V measured by an observer with four- velocity U^. We denote by e M and 
b v the electric and magnetic fields, respectively, measured by a comoving observer . 

2.2.1. Maxwell's equations. The equations governing the evolution of the electro- 
magnetic fields are Maxwell's equations, which can be written in terms of the Faraday 
tensor as *F^ V . v = and F^ v . v — AnJ^, where is the electric four-current. Under 
the assumption that Ohm's law is fulfilled, the latter reads J 11 = p q it M + ae^, where 
p q is the proper charge density and a is the electric conductivity. Maxwell's equations 
can be simplified if the fluid is a perfect conductor (a — > 00). In this case, to keep 
the current finite, e M must vanish. This case corresponds to the so-called ideal MHD 
condition. Under this condition the four-vector electric field E^ L can be expressed 
in terms of the four- vector magnetic field B^, and, thus, only equations for B l are 
needed. It is convenient to choose as observer the Eulerian observer, {7 M = n M , for 
which the temporal component of the electric field vanishes, E^ — (0, —SijkV^ B k ). In 
this case the first set of Maxwell's equations reduce to the divergence-free condition 
plus the induction equation for the magnetic field 




Fluid — 



phu^u v + Pg' 



TZi = ****** ~ -g^F x5 F xs , (4) 



ViB* 1 = 0, 



dB n 



(•5) 



dt 



Towards relativistic simulations of magneto-rotational core collapse 



4 



where B* % = yfjB 1 and v' n = av l = av l — j3 l , v l being the three- velocity as measured 
by the Eulerian observer. Here, 7 denotes the ratio of the determinants of the three- 
metric and the flat three- metric, 7 = 7/7. 



2.2.2. Magnetic flux conservation. The total magnetic flux through a closed surface 
A enclosing a volume V can be calculated as a surface integral of the "starred" 
magnetic field as 

$ T = / B* ■ dA = [ V ■ B*dV = 0, (6) 
JA=dV Jv 

after applying Gauss theorem and the magnetic field divergence-free constraint. 
Quantities in boldface correspond to three- vectors. The scalar product (•), and the 
cross product (x) used below, are defined with respect to the flat three- metric 7^. 
This equation implies that no source of magnetic flux exists inside the volume V and, 
therefore, the magnetic flux is a conserved quantity as = 0. 

In addition, if we consider a generic surface A (without the restriction of having 
to enclose a volume) , the time variation of the magnetic flux through the surface is 

where we have used the induction equation ([5]) and Stokes theorem to transform the 
surface integral into a line integral along the curve C enclosing A, and the equality 
E* = v* x B* . The two properties just inferred allow us to design a numerical 
algorithm to solve the induction equation and the divergence constraint in a way that 
ensures the conservation of the magnetic flux. 



dA = - j> E* ■ dl, (7) 



2.2.3. Conservation laws. In the ideal MHD limit the energy-momentum tensor of 
the electromagnetic field can be written in terms of the magnetic field measured 
by a comoving observer. The total energy-momentum tensor is thus given by 

(ph + b 2 ) u» ' 1 ■ ''' 



1 Fluid 



' EM 



P + — ) g» v - , (8) 



where b 2 — b^U L . We define the magnetic pressure P mag = 6 2 /2 and the specific 
magnetic energy e mag — b 2 /2p, as their effects on the dynamics are similar to those 
played by the pressure and the specific internal energy. The evolution of the magneto- 
fluid is determined by the conservation law of the energy-momentum T^.^ = and 
the continuity equation J^.^ = 0, for the rest-mass current J p = pu^. Following the 
procedure laid out in [TT], the magnetic field can be accounted for by choosing the 
conserved quantities in a similar way to the purely hydrodynamical case [15] 

D = pW, (9) 

S l = {ph + b 2 )W 2 v l - (10) 

/ 6 2- 
r = E- D = (ph + b 2 )W 2 -\P 



-)-a 2 (b°) 2 -D. 



(11) 



With this choice, the system of conservation equations for the fluid and the induction 
equation for the magnetic field can be cast as a first-order, flux-conservative, 
hyperbolic system, as 



1 



dt 



dx l 



Q, 



(12) 
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with the state vector, flux vector, and source vector given, respectively, by 

U = [D,Sj,T,B k ], (13) 



Dv\SjV l + 5} [P + 



bjB* 
2 J ' ~W 
b°B l 



Q = 



w 

d In a 



, v l B H 



v k B i 



'2 ' V ftc" 



r^J ,o' 



(14) 



(15) 



We note that these expressions contain components of the magnetic field measured by 
both, a comoving observer and an Eulerian observer. The two are related by 



WB l Vi 



B l + ab°u l 

w 



v = 



(16) 



2.2.4- Hyperbolic structure. Anticipating the numerical methods that we use to solve 
the conservation equations (112|) . we need to know the wave structure of the hyperbolic 
system of equations. The associated flux- vector Jacobians in every direction are 7x7 
matrices, and the solution of the eigenvalue problem [16j E] leads to seven types of 
waves, which may appear when solving the Riemann problem for each direction i: 
the entropic wave A* = av l — (3 l , the Alfven waves \ l a ± = (b l ± y/ ph + B 2 u l )/{b° ± 
\J ph + B 2 u°), and the magnetosonic waves. There is no analytic expression for the 
latter, which must be computed as the solution of a quartic equation (see [llj for 
details). Among the magnetosonic waves, the two solutions with maximum and 
minimum speeds are called fast magnetosonic waves A^ ± , and the two solutions in 
between are the so-called slow magnetosonic waves A* ± . The seven waves can be 
ordered as follows A}_ < A J Q _ < A*_ < A* < A* + < A^ + < \) + . 



2.2.5. The passive field approximation. In the collapse of weakly magnetized stellar 
cores, it is a good approximation to consider that the magnetic field entering in the 
energy-momentum tensor of © is negligible when compared with the fluid part, i.e. 
-Pmag "C P, e mag <C e, and that the components of the anisotropic term of T^ v satisfy 
b^b v <§C phu^u w + Pg^ v . With such simplifications the remaining system of equations 
comprises the hydrodynamics equations (with no magnetic field) and the induction 
equation. In this case, the magnetic field evolution does not affect the dynamics of 
the fluid, but it does depend on the evolution of the matter due to the presence of the 
velocity components in the induction equation. 

Such "test magnetic field" (or passive field) approximation is employed in the 
core collapse simulations reported in this paper. In this approximation the seven 
eigenvalues of the GRMHD Riemann problem reduce to three 

(17) 



\ l — \ l — \ l — \ l 

A hydro — A e — A a± — A s± ; 



A ± hydro ~ A /± ) 



where AQ hydro and Aj_ hydro are the eigenvalues of the Jacobian matrices of the 



hydrodynamics equations (see |17j). 



3. Numerical framework 



The GRMHD numerical code presented in this paper is based on the hydrodynamics 
code described in [TOj [4] , and on its extension discussed in [18] . The code performs 
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the coupled evolution of the equations governing the dynamics of the spacetime, fluid, 
and magnetic fields in general relativity. The equations are implemented in the code 
using spherical polar coordinates {t, r, 6, </?}, assuming axisymmetry with respect to 
the rotation axis and equatorial plane symmetry at = ir/2. For the various types 
of PDEs implemented in the code we use the most appropriate numerical methods, 
which will be described next. 

3.1. The hydrodynamics solver 

Since we adopt the passive magnetic field approximation, the fluid evolution is not 
affected by the magnetic field. Such evolution is performed using a high-resolution 
shock-capturing (HRSC) scheme which numerically integrates a subset of equations in 
system f|12[) . the one corresponding to the purely hydrodynamical variables (-D, Si,r). 
HRSC methods ensure numerical conservation of physically conserved quantities and a 
correct treatment of discontinuities such as shocks (see e.g. .17] and references therein). 
Following [11] we use a linear reconstruction procedure with a minmod slope limitcr, 
which yields second order accuracy in space. Correspondingly, the time update of the 
state vector U is done using the method of lines in combination with a second-order 
accurate Runge-Kutta scheme. The numerical fluxes at cell interfaces are obtained 
using "incomplete" approximate Riemann solvers, i.e. solvers that do not need the 
full characteristic information of the system. This kind of solvers are particularly 
useful in GRMHD, where the full set of eigenspeeds of the flux- vector Jacobians is not 
known in a closed form [TT] . We have implemented the HLL single-state solver of [T9] 
and the symmetric scheme of [20) (KT). Both solvers yield results with an accuracy 
comparable to complete Riemann solvers (with the full characteristic information), as 
shown in simulations involving purely hydrodynamical special relativistic flows [21] 
and general relativistic flows in dynamical spacetimes [52]. Test of both solvers in 
GRMHD have been reported recently by [TTj . 

3.2. Magnetic field evolution 

The magnetic field evolution needs to be performed in a different way to the rest of the 
conservation equations, since the physical meaning of the corresponding conservation 
equation is different. Although the induction equation can be written in a flux 
conservative way, a supplementary condition for the magnetic field has to be given 
(the divergence constraint), and it has to be fulfilled at each time iteration. As shown 
before, the physical meaning of these two equations is the conservation of the magnetic 
flux in a closed volume, in our case each numerical cell. Therefore, an appropriate 
numerical scheme has to be used which takes full profit of such conservation law. 
Among the numerical schemes that satisfy this property (see |23j for a review), the 
constrained transport (CT) scheme |24j has proved to be adequate to perform accurate 
simulations of magnetized flows. Our particular implementation of this scheme has 
been adapted to the spherical polar coordinates used in the code. 

3.2.1. CT scheme in spherical polar coordinates. In order to implement the CT 
scheme in our numerical code we have to analyze how the magnetic flux behaves at 
the surface of a numerical cell. Therefore, we apply ([7]) for the magnetic flux evolution 
in the interfaces along the direction of each coordinate (see figure [l}. To do this we 
assume that B* 1 is constant over each cell surface, and E* is constant along each cell 
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Figure 1. Schematic representa- 
tion of the numerical cell. The 
time derivative of the magnetic 
fluxes $i over the r interface, the 8 
interface, and the <p interface can 
be written as line integrals along 
the corresponding closed path a-b- 
c-d. Only the first case is shown in 
the plot. 




edge. In the r direction this yields 



dt 



- d 
AA r -B* r 
at 



[E;Ai v ] a - [e;aI 



(18) 



keeping in mind that the axisymmetry condition imposes that [Eg Alg]b — [EgAle]d- 
We define AAi = J dAi as the surface of the cell interface perpendicular to the i 
direction and AU = J dk as the length of the cell edge in the i direction. In the 9 and 
p directions similar relationships can be found. 

We represent all these quantities in the numerical grid, which has a total number 
of n r x ng points (cell centers) labelled (i j), with i = 1 . . . n r and j = 1 . . . n$. Cell 
interfaces between neighboring cells are denoted as (i + 1 / 2 j) for the radial ones and 
(i j + 1 / 2 ) for the angular ones. Hence, indices (i + 1 / 2 j + 1 / 2 ) denote cell edges in the p 
direction. In order to implement the CT scheme we do not define the poloidal magnetic 
field cell-centered but at the interfaces, i.e #*+y 2 j an d ^*/+i/ 2 - However, as we assume 
axisymmetry, we use a cell-centered toroidal field B**f because this component does 
not play any role in the CT scheme. Therefore, we consider hereafter the poloidal 
components of the magnetic field only. With such considerations the above equations 
for the magnetic flux time evolution lead to equations for the magnetic field evolution 
(CT scheme). The poloidal components read: 



dB* 



dt 



e;ai v 



i+72 J-72 



e; ai v 



*+72 i+72 



dB* 



dt 



»+72 3 



i 3 + 1 h 



e; ai v 



AA r 



j+ 1 / 2 



i+72 3 



e; ai v 



»-7a J+72 



AA 



i j+ 1 h 



(19) 



(20) 



The total magnetic flux through the cell interfaces is given by 

$T i j = $r i+7 2 j - $r i-i/ 2 j+®e i j+7 2 - *0 i j- 1 ^ ( 21 ) 

where we have taken into account that the total flux in the ip direction is zero owing 
to the axisymmetry condition. The time evolution of the total magnetic flux evolved 
with the CT scheme satisfies by construction that [9$ T /9f]j j = 0, and therefore 
every numerical scheme constructed on the basis of the CT scheme will conserve 
magnetic flux up to machine accuracy. If we are able to generate initial conditions 
which satisfy the divergence constraint, i.e. with <J>t = at each numerical cell, then, 
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such constraint will be preserved during the numerical evolution. The way to build 
such initial data is explained below. 

The next step is to choose a discretization of the integrals by making explicit the 
values of the interface surfaces and of the edge lengths 

r '+y 2 A(cos0) J Aip, (22) 

1 



Al r i j+ i /2 = Ar, 



AA r j+Va j = -rl 



Al 



8 i+Va 3 



AA 



e i 3+V2 



■ sin^- + i/ 2 Ar~ Aip, (23) 



Al 



r i+ i/ 2 sin ( 



where Ar, = 



- n-i/ 2 , Ae 3 = % +V2 - 



AA, 



5 Ar? Atf„ 



(24) 



"i-Vn = r i+y 2 - r ty 2 > a ( cos0 )j = 
cos0j + i/ 2 — cos#y_y 2 , and Ay is arbitrary due to axisymmetry (hence, the cell size on 
the ip direction does not play any role in the numerical scheme). Taking into account 
the above expressions the evolution equations for the poloidal magnetic field read 



dB* r 
dt 

dB* e 
dt 



j-V= E % 1+V2 3-V2 



»+V2 i 



V2 3 



A(cos( 



»+V 2 »+y 2 J+ y 2 r »-V2 % t -y 2 3 +y 
Ar? 



£7* 



= 2 



(25) 



(26) 



which are used in the numerical code to update the magnetic field. The only remaining 
aspect is to give an explicit expression for the value of E* . A practical way to calculate 
E* from the numerical fluxes in the adjacent interfaces |25j is 

% i+y 2 j+y 2 = ~4 ( Fr )i 3+V2 + ( pr )i+i i+y 2 _ ( F8 )I+y 2 j ~ (-01+ y 2 i+y > ( 27 ) 

where the fluxes are obtained in the usual way by solving Riemann problems at the 
interfaces. The combination of the CT scheme and this way of computing the electric 
field is called the flux-CT scheme. It is used in all numerical simulations reported in 
this paper. 



3.3. The metric solver 

The CFC metric equations are 5 nonlinear elliptic coupled Poisson-like equations which 
can be written in compact form as Au{x) = f(x;u(x)), where u — u k — (</>, a</>, j3 3 ), 
and / = f k is the vector of the respective sources. These five scalar equations for 
each component of u couple to each other via the source terms that in general depend 
on the various components of u. We use a fix-point iteration scheme in combination 
with a linear Poisson solver to solve these equations. Further details on such metric 
solver can be found in [HI E] • 



4. Initial magnetic field configurations 

Since the CT scheme only preserves the value of V • B* but does not impose the 
divergence constraint of the magnetic field itself, we have build initial conditions 
which satisfy such condition. To do this we calculate the initial magnetic field from 
a vector potential A*, such that B* — V x A*. We choose the following numerical 
discretization of this equation 

i sitlQa . i/„ A* . , ,, ■ , — sin#.,_ i/„ A* . , ,, - ,, 

r>*r _ _l 3+ / 2 ip 1+1/2 j+yja j g <p '+y 2 j-y 2 / 08 ^ 

B< +V2i- — A(cos^ ' (28) 
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Figure 2. Global error 
in the toroidal magnetic 
field, <t, after a time evolu- 
tion of 1 ms for tests TTA 
(x), TTB (O ), and TTC 
(+) as a function of 1/f 
for a sequence of models 
with grid resolutions 80 X 
10 (/ = 1), 160 x 20 (/ = 
2) and 320 x 40 (/ = 4). 
The solid lines are the best 
fits of each test to a power 
law. 



0.0001 



le-06 




u *e r »+V2 A * v »+y 2 j+y 2 r »-V 2 A % t -y 2 3 +y 2 , . 

B i 3+V2 = ~ 2 ■ (29) 

which warrants that the numerical magnetic flux $x over a cell is zero (up to round-off 
error) in the specific numerical representation of our CT scheme. 

For our code tests and core collapse simulations we use two possible poloidal 
magnetic field configurations as initial conditions: a) The homogeneous "starred" 
magnetic field, in which B* is constant and parallel to the symmetry axis, and b) the 
magnetic field generated by a circular current loop of radius r mag [26] . Note that in 
both cases, we employ the "starred" magnetic field, since the divergence constraint is 
valid for this quantity when computed with respect to the flat divergence operator. In 
this way we can extend any analytic prescription for the magnetic field given in flat 
spacetime in an easy way. Note that in the presence of strong gravitational fields the 
magnetic field B is deformed with respect to B* due to the curvature of the spacetime, 
although the divergence constraint is automatically fulfilled. 



5. Code tests 



We have designed several tests in order to check the accuracy of our numerical code 
when solving the induction equation with the numerical methods described in the 
previous sections. The "toroidal test" is designed to assess the ability of the code to 
maintain equilibrium magnetic field configurations (labelled TTA and TTB) and to 
correctly compute the amplification of the toroidal magnetic field as it is wound up 
by a rotating fluid (TTC). On the other hand, the "poloidal test" (PT) is designed 
to check whether the code can correctly compute the compression of the poloidal 
magnetic field in a spherical collapse, and its ability to handle the presence of radial 
shocks. 



5.1. Toroidal test 

We consider a rotating stationary configuration with no meridional flows, v* r = v* e = 
0, and v* v = fi*(r, 0) r sin#, where f2*(r, 6) stands for the rotation law. Under 
these conditions and in the passive field approximation, the induction equation can 
be integrated analytically. The solution shows that the poloidal component of the 
magnetic field remains constant and the toroidal component grows linearly with time 
as B* v (t) = B* v (t = 0) + 1 rsinOB* ■ Vft*. In order to test whether the numerical 
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Figure 3. Local order of convergence (color coded) for the TTA test after a 
total time evolution of 1 ms. White color is used for values larger or equal than 
3.0. The axes represent the number of cells of the reference grid in the radial and 
angular direction. 



code is able to recover this solution, we investigate three particular cases. For all 
three tests we consider a non-evolving fluid of constant density p — 10 14 g cm -3 
filling the computational domain, and an equally-spaced grid in both the angular and 
radial directions. The outer radial boundary is located at 20 km. We also assume 
a static background spacetime with a flat metric. We note that our choice of units 
is arbitrary since the background metric is flat. However, we maintain the adopted 
units to keep a link with the scenario for which the code is prepared, namely core 
collapse. Test TTA consists in a rigidly rotating fluid, Q.* — fi*, and the magnetic 
field generated by a circular current loop of radius r mag = 6 km. Test TTB includes a 
homogeneous magnetic field and a differentially rotating fluid with a rotation law given 
by n*(r,0) = (A 2 tt*)/(A 2 + (rsinfl) 2 ), with A = 6 km. In both cases the analytic 
solution corresponds to a non-evolving toroidal magnetic field because B* ■ VfT = 0. 
The third test, TTC, consists in a homogeneous magnetic field and a fluid satisfying 
a rotation law depending on the radius f2*(r) = (A 2 fl*)/ (A 2 + r 2 ), and A = 6 km. In 
this case the analytical solution is given by 

= + S[U29 - (30) 

In all three cases the central angular velocity fl* is such that the center rotates 10 
times during the complete evolution, namely 1 ms. Once this time is reached we 
compute the local error cry = \B V — Bf t ^. \ between the analytical solution and 
the numerical one, and the global error a calculated as the L2-norm of the local error. 

If we consider a reference (low resolution) grid n r rc f x ng re f with an error cr re f and 
a subsequent grid with a factor / higher resolution, i.e. n r = fn rTe f and ng = fng Te {, 
then the error on the new grid, <r/, is related to the error on the reference grid as 
o~f = <j re f (l/f) N , where N is the order of the numerical scheme. By fitting the values 
of cr/ as a function of l/f we can easily calculate the convergence order N. For our 
test simulations we consider a reference grid of 80 x 10 cells (/ = 1) , and higher 
resolution grids of 160 x 20 (/ = 2) and 320 x 40 (/ = 4) cells, respectively. Figure [2] 
shows the resulting a/ for the three tests versus l/f. Our results show that (i) the 
order for the TTC test (N — 1.54, obtained by fitting the data to a power law) is 
smaller than for the TTA and TTB tests (N = 2.35 and 2.64, respectively), and (ii) 
the order for cases TTA and TTB is TV > 2, and hence higher than the theoretical 
expectation (second order given by the time and space discretization order). 

The main difference between the case TTC and the cases TTA and TTB, which 
explain these results, is that in the first case there is a component of the magnetic 
field, B* v , which grows linearly in time, while in the other two cases no components 




Figure 4. Results for the poloidal test (PT) for different grid resolutions. The 
left panel shows the global error as a function of time during the infall phase, 
while the right panel shows the quantity rD*/B versus m at t = 35 ms. The 
location of the inner core of the PNS (dashed vertical line) and the location of 
the shock (dotted vertical line) for the highest resolution model are also shown. 



evolve. Hence, the order of convergence for tests TTA and TTB is higher than for test 
TTC. We suspect that this is due to more precise numerical cancellations in Eqs. (flZij]) 
and (gU) (which have E* ^ 0) in tests TTA and TTB. This can be explained by 
computing the local order of convergence, i.e. the order obtained when computing the 
errors of each numerical cell, er^ , instead of the global error a. The results for test 
TTA are displayed in figure [3] (similar plots can be obtained for the other two cases). 
It can be seen that at some particular grid zones the order of convergence is larger 
than two, while at most locations it remains around two. 

We note that the computation of the initial magnetic field for test TTA is not 
analytic but the result of an expansion with infinite terms. This expansion is computed 
to a required level of accuracy, typically smaller than the accuracy required during the 
evolution in the routine to recover the primitive variables from the state vector. For 
a particular direction, however, close to the center of the "peculiar" cells in figure [31 
the solution needs infinite terms, and the resulting solution is thus less accurate at 
these points. The physical meaning of this direction in the equatorial plane is that 
it coincides with the radius of the current loop which generates the initial magnetic 
field. 



5.2. Poloidal test 

Next, we consider a test in which only radial velocities of the fluid are allowed, i.e. 
v* r 7^ and v* e = v* v = 0. We also consider initially a purely poloidal magnetic field. 
In this case, it can be easily shown from the induction equation ([5]) and the continuity 
equation (fT2)) that the following equivalence holds in the equatorial plane 

Ft (^) = ITr ( r !^) ' (31) 

which is an advection equation for rD*/B* e with velocity v* r at the equator. Since 
only radial velocities are allowed, it is possible to define a Lagrangian coordinate 
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Figure 5. Time evo- 
lution of the magnetic 
energy parameters /3 m ag 
(solid line), j3 v (dashed 
line), and /3 po i (dot- 
dashed line) during 
magneto-rotational col- 
lapse. The vertical dotted 
line indicates the time of 
bounce. 
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system in which the value rD*/B does not change with time. The easiest way of 
checking this condition is to calculate [r D* / B 6 *] (m), i.e. as a function of the mass 
enclosed within a radius r defined as m(r) = 4n r' 2 dr' D* (r'). This quantity should 
remain constant in time, provided the magnetic field is correctly evolved with our 
numerical code. 

To test this fact we perform simulations of spherical collapse in which the above 
conditions are satisfied. The initial model is a spherical equilibrium 4/3-polytrope 
(TOV) with central density p c — 10 10 g cm" 3 and a homogeneous (starred) magnetic 
field. We use a hybrid EOS (see [10]). We induce the collapse by reducing the initial 
adiabatic index to 71 = 1.28. As the EOS stiffens at nuclear matter density the 
star bounces (tb = 30 ms) and a shock forms, which travels outward. We calculate 
the global error a in the collapse phase as the L2-norm applied to the differences 
in the quantity [r D* j B e *~\ (m) between its initial value and its value at subsequent 
time steps. In the left panel of figure 2] we show the evolution of such error during the 
collapse for different (r, #)-grid resolutions, equally-spaced in the angular direction and 
logarithmically spaced in the radial direction, except for the inner 20 radial grid points 
which are equally-spaced. In each case the errors are below 1%, even for the coarsest 
grid, and the computed order of convergence during the collapse (at t = 20 ms) is 
1.41. 

In order to check our numerical code after core bounce we plot in the right panel 
of figure [1 [rD*/B 6 *] (m) at time t = 35 ms, for the three different grid resolutions. 
The initial profile is shown by the dotted curve. The effect of the travelling shock is 
seen as a small spike, that becomes narrower as the resolution is increased. On the 
other hand, deviations from the initial profile occur near the inner core boundary of 
the PNS (where the density exceeds nuclear matter density). This numerical error 
becomes smaller as the resolution covering the inner region is increased by using a 
grid of 400 x 60 zones, with 100 equally-spaced radial zones in the inner part instead 
of 20. This resolution is the same as the one used in magneto-rotational core collapse 
simulations. It guarantees sufficient resolution in the region where the PNS forms. 
We have performed comparisons of the HLL approximate Riemann solver and the KT 
symmetric scheme, finding almost identical results (in agreement with [211 [22l [TTj V 
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Figure 6. Configuration of the innermost region of the star at the end of the 
magneto-rotational collapse of its core (t = 60 ms). The panel shows the logarithm 
of density, logp, in g cm -3 , the distribution of the velocity field, v l (arrows), and 
contours of the specific internal energy. All axes are in km. 



6. Magneto-rotational core collapse 

We end the paper by discussing a representative magneto-rotational core collapse 
simulation. A parameter-space survey of a large sample of models will be presented 
elsewhere. The initial model, which we label A1B3G3-D3M0, is a general relativistic 
and rotating 4/3-polytrope in equilibrium, and is the magnetized version of the purely 
hydrodynamical model A1B3G3 of [10], which is differentially and rapidly rotating. 
The configuration of the initial magnetic field is a circular current loop of radius 
r mag = 400 km with a strength at the center of the loop of B c = 10 10 G. For such value 
the passive field approximation is adequate, since, locally, P mag /P < 10~ 6 throughout 
the simulation. The collapse begins by lowering the adiabatic index to ji = 1.31. The 
same hybrid EOS as in the PT test is used. We employ the numerical grid described 
in the previous section. 

Figure [5] shows the evolution of the energy parameter for the magnetic field, 
/3ma g , defined as the ratio of magnetic energy, E mag = | J d 3 xWb 2 , to the potential 
energy E pot (see e.g. [ID])- In order to analyze the growth of the magnetic field, we 
separate the effect of the different components of the magnetic field into (3 V , for the 
toroidal component, and /3 po io = Anag — Pip, for the poloidal component, which are 
also plotted in the figure. As the collapse proceeds, and leaving aside the effects of 
the magneto-rotational instability (MRI), the magnetic field grows by at least two 
reasons: First, the radial flow compresses the magnetic field lines, amplifying the 
existing poloidal and toroidal magnetic field components. Second, during the collapse 
of a rotating star differential rotation is produced, even for rigidly rotating initial 
models, as shown in numerical simulations (e.g. [101 [S]). This fact can be easily 
explained considering that energy and angular momentum are roughly conserved in 
cylindrical regions (see e.g. [27]). Hence, if a seed poloidal field exists, the fi-dynamo 
mechanism acts winding up poloidal field lines into the toroidal component. This 
(linear) amplification process generates a toroidal magnetic field component, even 
from purely poloidal initial configurations. The toroidal component of the magnetic 
field is affected by the two effects while the poloidal field is only amplified by the first 
effect (radial compression). Thus, even if the initial magnetic field configuration is 
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Figure 7. Configuration of the innermost region of the star at the end of the 
magneto-rotational collapse of its core (t = 60 ms). The left panel shows the 
logarithm of the poloidal component of the magnetic field, log |-B po i |, in G and 
the magnetic field lines in the r-8 plane (lines). The right panel shows, B^/[B po i |. 
All axes are in km. 



purely poloidal, the toroidal component dominates after some dynamical time. This 
can be seen in figure [5] as {3$ (dashed line) grows much faster than /3 po io, particularly 
after bounce (t ~ 50 ms) when the radial compression mechanism stops. We note that 
as the magnetic field considered is weak enough not to affect the dynamics, the final 
/?mag is much less than unity. 

Figures [6] and [7] show the innermost 100 km of the star after core bounce. A 
PNS has formed at the center with two distinct regions. The inner core of ~ 10 
km, where nuclear matter density is reached, has a mixed poloidal-toroidal magnetic 
configuration, while there exists a surrounding shell extending to ~ 50 km, with sub- 
nuclear densities, that has a dominant toroidal magnetic field which grows linearly due 
to the J7-dynamo mechanism resulting from the strong differential rotation of this shell. 
From the linear growth of this component we estimate that the magnetic field is likely 
to reach saturation values of ~ 10 15 G on a timescale of several seconds. However, 
by estimating the timescale of the fastest growing unstable mode of the MRI in our 
simulation, we have checked that a significant fraction of the newly-formed PNS, as 
well as the region behind the shock at the moment of its formation, could be affected 
by the MRI. Such instability can grow on dynamical timescales of several tens of ms 
[9] , much faster than the f2-dynamo mechanism. However, it is still an open issue |28j 
whether, for the weak magnetic fields encountered in core collapse progenitors, the 
amplification of the magnetic field is going to become important for the dynamics. If 
that were the case, the passive field approximation would not be valid in such regions, 
becoming necessary a full magnetic field treatment. 

7. Summary 

We have presented a new general relativistic hydrodynamics code specifically designed 
to study magneto-rotational, relativistic, stellar core collapse. The code is built on 
an existing hydrodynamics code which has been thoroughly applied in the recent 
past to study relativistic rotational core collapse [4]. It is based on the conformally- 
flat approximation of Einstein's field equations and conservative formulations for the 
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magneto-hydrodynamics equations. As a first step towards magneto-rotational core 
collapse simulations, the code assumes a passive (test) magnetic field, a justifiable 
assumption since weakly magnetized fluids are present in this scenario . The paper 
has focused on the description of the technical details of the numerical implementation, 
with emphasis on the magnetic field module. A number of tests have been presented 
and discussed, and the convergence properties of the code have been analyzed. 
In addition, a representative magneto-rotational core collapse simulation has been 
presented. The amplification of the magnetic field due to the fi— dynamo mechanism 
has been estimated (acting on timescales of several seconds). More detailed studies 
of faster amplification mechanisms (namely, MRI) are necessary for initially weak 
magnetic fields, since they could dominate the post-bounce dynamics within a few ms 
and might have a major effect on the gravitational wave signal. 
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